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Abstract  Intrinsic  neuronal  and  circuit  properties  control  the  responses  of  large  ensembles 
of  neurons  by  creating  spatiotemporal  patterns  of  activity  that  are  used  for  sensory  process¬ 
ing,  memory  formation,  and  other  cognitive  tasks.  The  modeling  of  such  systems  requires 
computationally  efficient  single-neuron  models  capable  of  displaying  realistic  response 
properties.  We  developed  a  set  of  reduced  models  based  on  difference  equations  (map-based 
models)  to  simulate  the  intrinsic  dynamics  of  biological  neurons.  These  phenomenological 
models  were  designed  to  capture  the  main  response  properties  of  specific  types  of  neurons 
while  ensuring  realistic  model  behavior  across  a  sufficient  dynamic  range  of  inputs.  This 
approach  allows  for  fast  simulations  and  efficient  parameter  space  analysis  of  networks 
containing  hundreds  of  thousands  of  neurons  of  different  types  using  a  conventional 
workstation.  Drawing  on  results  obtained  using  large-scale  networks  of  map-based  neurons, 
we  discuss  spatiotemporal  cortical  network  dynamics  as  a  function  of  parameters  that  affect 
synaptic  interactions  and  intrinsic  states  of  the  neurons. 
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1  Introduction 

The  most  realistic  approach  for  simulation  of  neuronal  behavior  is  based  on  Hodgkin- 
Huxley-type  models  [1]  where  each  ionic  current  is  described  using  voltage-dependent 
opening  and  closing  rates  for  the  gating  variables.  Examples  of  successful  Hodgkin- 
Huxley-type  models  include  cortical  cells  [2-6],  thalamic  relays  [7-9],  reticular  cells 
[10-12],  hippocampal  neurons  [13,  14],  and  antennal  lobe  neurons  [15,  16],  among  others. 
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Unfortunately,  the  high  dimensionality  and  complexity  of  nonlinear  functions  that  constitute 
Hodgkin-Huxley-type  models  hamper  their  application  to  the  simulation  of  large-scale 
networks  containing  a  realistic  number  of  neurons. 

Integrate  and  fire  models  (IF)  [17-19]  are  a  class  of  models  commonly  used  in  large- 
scale  simulations.  In  the  IF  model,  a  neuron  simply  integrates  its  inputs  and  generates  a 
spike  when  a  threshold  is  reached.  After  the  spike,  the  membrane  voltage  is  reset  to  zero, 
followed  by  a  refractory  period  during  which  spiking  is  impossible.  The  disadvantage  of 
IF  models  is  that  the  firing  patterns  are  oversimplified  and  do  not  describe  the  range  of 
experimental  data.  IF  models  also  fail  to  reproduce  nonlinear  intrinsic  resonance  properties 
of  biological  neurons.  Generalizations  of  this  simple  model  include  leaky  IF  models  that 
introduce  a  leak  term  in  the  dynamics  of  the  subthreshold  membrane  voltage  and  the 
integrate-and-fire-or-burst  neuron  model  that  reproduces  rebound  bursting  observed  in  some 
cell  types  [20,  21].  However,  generalizations  of  the  IF  model  that  help  describe  realistic 
firing  patterns  significantly  increase  the  complexity  of  the  IF -based  model  and  reduce 
the  speed  of  computer  simulations.  One  recent  exception  proposed  by  Izhikevich  [22]  is 
designed  in  the  form  of  a  two-dimensional  system  of  ordinary  differential  equations  that 
replicates  a  variety  of  physiological  firing  patterns.  For  computer  simulations,  this  model 
was  rewritten  as  a  system  of  difference  equations  using  Euler’s  integration  scheme  with  a 
step  size  small  enough  to  allow  the  resolution  of  spikes. 

We  have  proposed  a  way  to  improve  the  computational  efficiency  of  network  simulations 
by  designing  a  model  of  spiking  neurons  that  neglects  the  duration  and  shape  of  individual 
spikes  but  captures  the  dynamics  occurring  over  other  time  scales  [23-25].  This  is  achieved 
by  designing  the  dynamical  system  in  the  form  of  difference  equations  (i.e.,  a  map).  The 
explicit  discretization  in  time  permits  a  relatively  low  sampling  rate  of  0.5  ms,  enabling 
efficient  simulation.  However,  in  contrast  to  continuous  time  models  simulated  at  a  low 
sampling  rate,  each  spike  is  given  by  a  single  sample  in  the  waveform  data  and  is,  therefore, 
never  missed.  In  addition  to  the  firing  patten,  other  relatively  slow  dynamical  processes 
are  captured  correctly  by  the  map-based  model  because  they  are  unaffected  by  the  temporal 
discreteness  of  the  model.  Attractive  features  of  this  map-based  model  include  the  simplicity 
of  the  equations,  the  ability  to  describe  a  broad  range  of  firing  patterns  found  in  biological 
neurons,  and  the  opportunity  to  adopt  models  of  synaptic  inputs  used  in  simulations  of 
networks  with  Hodgkin-Huxley  equations. 

In  this  paper,  we  start  by  introducing  the  map-based  model.  The  dynamics  is  discussed 
as  a  function  of  the  model  parameters.  In  the  second  part  of  the  manuscript,  we  use  this 
generic  approach  to  design  specific  models  describing  different  classes  of  thalamic  and 
cortical  neurons  that  are  analyzed  in  response  to  constant  and  periodic  external  inputs.  In 
the  last  section,  we  apply  this  approach  to  study  intrinsic  mechanisms  and  synchronization 
properties  of  fast  gamma-range  (30-80  Hz)  oscillations  in  large-scale,  two-dimensional 
cortical  network  models. 


2  Map-based  Approach  to  Modeling  a  Typical  Action  Potential 

An  action  potential  (spike)  represents  one  of  the  best  known  characteristic  features  of 
biological  neurons  and  plays  an  important  role  in  information  transfer  across  neuronal 
systems.  In  many  classes  of  biological  neurons,  action  potentials  are  generated  by  the 
dynamics  of  voltage-dependent  Na+  and  K+  conductances  [1].  In  this  section,  we  will 
introduce  a  reduced  model  of  a  biological  neuron  that  simulates  action  potential  generation. 
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The  waveform  of  the  action  potential  can  be  replicated  using  the  reduced  model  written  in 
the  form  of  a  one-dimensional  nonlinear  difference  equation  (map)  as  follows: 


%n+ 1  —  folium  ^n)  •  (1) 

Here,  the  dynamical  variable  xn  represents  the  membrane  potential  of  a  biological  neuron. 
It  samples  the  action  potential  waveform  at  the  discrete  moments  of  time  n.1  The  time- 
dependent  variable  un  =  y+  includes  two  components:  the  parameter  y  defines  the 
baseline  (resting  potential)  of  the  neuron  and  the  time-dependent  term,  f3n  (index  n  indicates 
value  of  /3  at  time  n),  represents  the  external  input  to  the  model  (e.g.,  external  and  synaptic 
currents). 

The  shape  of  each  spike  of  neuronal  activity  is  replicated  here  using  a  piece-wise 
nonlinear  function  fa(xn,  un)  suggested  in  [24]: 


fa  (p^ri  ?  Mn  ) 


Qf/(1  —  Xn)  +  Un,  Xn  <  0 

a  +  un,  0  <  xn  <  a  +  un  and  xn-\  <  0 

—  1,  xn  >  a  +  un  or  xn-\  >  0. 


(2) 


The  control  parameter  a  in  (1)  is  usually  selected  from  the  range  3  <  a  <  4  for  tonically 
spiking  (nonbursting)  neurons  (see  more  about  this  below).  The  nonlinear  interval  of  the 
function  at  negative  values  of  xn  replicates  the  subthreshold  state  of  a  neuron  (e.g.,  the 
phase  of  the  resting  potential,  stimulation,  and  the  rising  phase  of  a  spike).  The  other  two 
conditions  of  (2)  are  involved  in  shaping  the  tip  of  the  spike  (made  by  a  single  sample  at 
xn  =  a  +  un)  followed  by  repolarization,  which  is  done  here  with  one  iteration  of  the  map. 
In  analogy  with  conductance-based  models,  these  two  intervals  represent  the  activation  of 
the  voltage-dependent  Na+  and  K+  conductances.  Note  that,  though  the  main  part  of  the 
map  (i.e.,  the  interval  xn  <  0)  is  a  one-dimensional  map,  the  condition  for  repolarization 
requires  the  use  of  xn-\ .  These  values  ensure  the  stability  of  the  repolarization  phase  in  the 
presence  of  external  stimuli  applied  to  the  model  (see  [24]  for  details). 

The  model  (1)  exhibits  two  main  dynamic  states  of  a  biological  neuron,  a  silent  state  and 
a  periodic  sequence  of  action  potentials.  Indeed,  depending  on  the  values  of  y  (assuming 
that  =  0)  (1),  the  map  may  have  a  pair  of  fixed  points  xps,  xpu  and  a  limit  cycle  (see  [24] 
for  details).  The  stable  fixed  point  xps  =  F(a,  un)  represents  the  resting  state,  and  the 
limit  cycle,  a  spike  sequence.  Remember  that  the  variable  un  includes  a  time- dependent 
component  fin  representing  input  to  the  system.  When  a  stimulus  (e.g.,  an  external  or 
synaptic  current)  is  applied  to  (1),  becomes  zero.  If  >  0  becomes  large  enough,  the  xps 
and  xpu  merge  and  disappear,  triggering  the  rising  phase  of  a  spike.  The  value  of  the  variable 
xn  where  the  fixed  points  merge  is  given  by  xth  =  1  —  «Joi.  The  corresponding  value  of  un 
is  given  by  u±  =  1  —  2^/a.  Thus,  like  conductance-based  models,  neurons  that  are  more 
depolarized  at  rest  (large  values  of  y)  require  less  external  input  (/ % )  to  reach  the  threshold 
for  spike  generation. 

Since  (1)  is  written  in  a  dimensionless  form,  it  generates  waveforms  with  xn  values 
ranging  from  —2  to  2.  To  recalibrate  the  waveforms  of  the  map  to  millivolts,  allowing 


1  More  often  than  not,  neuronal  dynamics  is  simulated  using  differential  equations  rather  than  difference  equa¬ 
tions  (for  example,  the  Hodgkin-Huxley  model  [1],  the  IF  model  [17-19],  and  the  Izhikevich  model  [22]). 
However,  it  is  important  to  emphasize  that  to  solve  a  differential  equation  numerically,  it  is  usually  rewritten 
in  the  form  of  a  difference  equation  (e.g.,  using  the  Euler  scheme,  the  time  derivative  dV/dt  is  replaced  by 
(Vn+i  -  V„)/r). 
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a  comparison  with  Hodgkin-Huxley  models  or  with  experimental  data,  one  can  use  the 
following  relationship: 


-50  mV 

vn  [mV]  =  - - —  X  Xn, 

1  -  v<* 


(3) 


where  the  denominator  stands  for  the  triggering  state  that  typically  occurs  at  a  depolarization 
level  near  —50  mV  This  normalization  sets  the  peak  of  the  spike  that  occurs  at  xn  =  a  + 
«th  =  (1  —  Va)2,  to  50  mV(^a  —  1),  which  reaches  50  mV  for  the  case  a  =  4. 

Thus  far,  we  discussed  the  model  with  no  explicit  time  scale.  Again,  using  an  analogy 
with  Euler’s  method  of  approximating  a  differential  equation,  the  time  interval  between  two 
successful  iterations  of  the  map  (1)  is  set  to  be  equal  to  the  numerical  integration  time  step. 
However,  since  the  model  (1)  is  introduced  as  a  difference  equation,  we  need  to  specify  the 
time  scale.  This  can  be  done  using  width  of  a  sodium  spike  that  is  constrained  by  the  kinetics 
of  Na+  and  K+  conductances  and  typically  falls  between  1  and  2  ms  (at  half  amplitude). 
Based  on  this  information,  we  can  set  the  time  interval  between  iterations  of  a  map  (1) 
as  0.5  ms. 


3  Slow  Dynamics  in  the  Map-based  Model  of  Neurons 

The  model  (1)  introduced  in  the  previous  section  is  capable  of  a  transition  between  silence 
and  tonic  spiking.  To  introduce  more  complex  dynamics,  including  spike  adaptation  or 
bursting,  we  need  to  introduce  a  slow  dynamical  variable,  which  may  also  be  described 
by  a  difference  equation.  In  this  section,  we  consider  an  approach  that  replicates  the  slow 
oscillatory  properties  of  a  neuron.  This  approach,  suggested  in  [24],  consists  of  a  map-based 
model  that  is  capable  of  generating  various  types  of  spiking  and  spiking-bursting  activity. 
The  system  including  additional  slow  dynamics  may  be  written  as  follows: 

%n+\  —  folium  "h  fin), 

yn+ 1  =  yn  -  +  1)  +  tier  +  non,  (4) 

where  xn  is  the  fast  and  yn  is  the  slow  dynamical  variable.  The  slow  time  evolution  of  yn 
is  achieved  by  using  small  values  of  the  parameter  fi  <<  1 .  Here,  the  control  parameter  y 
used  in  (1)  as  the  last  argument  of  the  function  fa  is  substituted  with  a  slow  variable  yn. 
Input  variables  f$n  and  on  introduce  the  action  of  synaptic  current  7syn  and  other  external 
currents  injected  into  the  neuron.  The  parameter  o  defines  the  resting  potential  of  the  model 
neuron. 

Detailed  analysis  of  the  individual  dynamics  of  model  (4),  with  =  0  and  crn  =  0, 
as  a  function  of  the  control  parameters  o  and  a ,  was  done  in  [24,  26].  In  the  two- 
dimensional  map-based  model,  the  single  fixed  point  O o  defines  the  resting  potential  level, 
which  is  now  set  using  the  control  parameter  o  instead  of  y,  used  in  (1).  When  a  < 
(7th  =  2  —  vW(l  —  />0>  ^e  system  converges  to  a  stable  fixed  point  Oo  with  coordinates 
xo  =  —  1  +  g  and  yo  =  —1  +  a  —  a/(2  —  a),  and  the  neuron  stays  silent.  At  the  threshold 
a  —  (7th,  the  fixed  point  loses  stability  due  to  a  subcritical  Neimark-S acker  bifurcation  (see, 
for  example,  [27])  and,  for  a  >  ath,  the  map  generates  spikes.2  Note  that  the  value  of  xn 


2  One  can  simplify  the  condition  for  the  threshold  by  introducing  a  new  control  parameter  crnew  =  o  —  2  + 
VW(1  —  l1)-  In  this  case,  the  threshold  value  occurs  at  7t^ew  =  0.  In  this  paper,  we  do  not  use  this  change  of 
the  control  parameter. 
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corresponding  to  the  threshold  level  fixed  point  is  xth  =  1  —  V<*(  1  —  l1)  ^  1  —  \[ol,  which 
is  about  the  same  as  in  the  model  (1)  and,  therefore,  normalization  (3)  remains  valid  for  this 
case. 

We  would  like  to  emphasize  that  the  addition  of  the  slow  variable  yn  to  the  model  (1) 
resulted  in  a  change  in  the  type  of  neuronal  excitability.  The  case  (1)  describes  a  type-1 
neuron  where  the  transition  from  silence  to  spiking  occurs  via  a  saddle-node  bifurcation. 
The  map  (4)  describes  a  type-2  neuron  where  the  transition  to  spiking  occurs  via  an 
Andronov-Hopf  bifurcation.  The  details  of  this  classification  may  be  found  elsewhere  [28]. 

The  bifurcation  diagram  of  the  map  (4)  plotted  in  the  parameter  plane  (a,  a)  is  presented 
in  Fig.  1.  The  diagram  shows  the  parameter  regions  corresponding  to  three  regimes  of 
behavior:  silence,  continuous  spiking,  and  generation  of  the  bursts  of  spikes.  Near  the 
border  between  the  spiking  and  spiking-bursting  activity,  shown  by  curve  Zts,  there  exists 
a  parameter  domain,  outlined  by  dashed  lines,  where  chaotic  regimes  of  continuous  spiking 
and  the  bursting  activity  occur;  see  [26]  for  details. 

The  diagram  also  shows  that  the  parameter  o  simulates  the  effect  of  a  depolarizing 
current  injected  into  the  real  neuron  [29].  To  get  a  better  view  of  this  similarity,  consider 
the  reaction  of  the  map  to  a  slowly  increasing  value  of  a.  The  test  with  a  slow  ramp  of 
injected  current  is  frequently  used  in  biological  experiments.  The  plots  of  the  map  reaction 
obtained  in  this  test  with  a  slew  rate  A o  =  0.00005  per  iteration  are  presented  in  Fig.  2. 
The  two  top  panels  of  this  figure  illustrate  a  transition  from  silence  to  spiking  through 
the  phase  of  spiking  bursting  activity.  The  transition  from  silence  to  continuous  spiking 
without  bursting,  which  takes  place  for  a  <  4.0,  is  shown  in  Fig.  2c.  Qualitatively  different 
transitions  (compare  Fig.  2c  with  Fig.  2a  and  b)  produced  by  the  model  are  observed  in  real 
neurons;  see,  for  example,  the  responses  of  the  PD  and  GM  neurons  of  the  stomatogastric 
ganglion  of  the  Californian  lobster  to  a  depolarizing  ramp  current  [30]. 

Notice  that,  due  to  a  continuous  increase  in  the  value  of  parameter  a,  the  oscillations 
in  the  map  occur  later  then  one  would  expect  based  on  the  threshold  values  ath  given  in 


Fig.  1  Diagram  of  dynamical  g  q 

regimes  produced  by  map  (4). 

Irregular  (chaotic)  spiking  and 
spiking-bursting  oscillations  are 
typical  within  the  intermediate 
region  outlined  by  dashed  lines 
and  crth 
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Fig.  2  Restructuring  of 
oscillation  in  map  (4)  with 
a  a=5.0,  b  a=4.5,  and  c  a=4.0, 
as  the  value  of  o  increases  with 
the  constant  slew  rate 
A o  =  0.00005  per  iteration 
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Fig.  1.  Usually,  spiking  in  biological  neurons  occurs  more  readily  for  fast  increases  in  the 
depolarization  current  than  for  the  slow  ones.  A  method  to  obtain  realistic  responses  from 
the  map  (4)  to  the  dynamically  changing  input  currents  is  discussed  in  the  next  section, 
where  we  will  describe  the  effects  of  injecting  a  current  In  using  input  variables  /3n  and  on. 
Note  that  on  in  (4)  adds  to  parameter  o  and,  therefore,  acts  similarly. 

3.1  Modeling  the  Response  to  Input  Currents 

Inputs  to  the  model  (4)  are  described  by  two  variables,  /3n  and  on ,  that  are  functions  of 
external  (e.g.,  synaptic  or  any  other)  currents.  an  acts  through  the  slow  subsystem  of  (4).  It 
changes  the  location  of  the  fixed  point  O  and  the  system  responds  to  it  by  slowly  moving 
towards  a  new  state.  Changes  of  the  value  of  /3n ,  on  the  contrary,  act  instantly  through  the 
fast  subsystem  of  (4)  and  quickly  change  the  state  of  the  neuron.  However,  after  the  changes 
in  /3n ,  the  system  (4)  slowly  drifts  towards  the  original  dynamics.  Using  the  change  of  the 
slow  variable  in  (4),  one  can  combine  both  these  input  parameters  and  rewrite  (4)  in  the 
form 


%n+ 1  —  Un), 

un+ 1  =  Un  -  /l(xn  +  1)  +  /Z<7  +  T)n, 

where  un  =  yn  +  /3n  is  the  new  variable  and  rjn  =  /icrn  —  /3n+\  +  /3n  is  the  new  input 
parameter.  Here,  one  can  see  that  the  map  reacts  only  to  the  changes  (derivative)  in  /3n, 
while  the  absolute  value  of  the  /3n  does  not  affect  the  transient  dynamics  of  the  map. 

In  practice,  it  is  convenient  to  use  the  system  in  the  form  (4).  Both  variables  /3n  and  an 
are  useful  in  modeling  a  variety  of  response  dynamics.  In  constructing  a  model  that  mimics 
the  response  of  a  real  biological  neuron,  one  should  experiment  with  the  map  and  define  a 
proper  balance  between  these  two  functions  of  external  current  to  achieve  the  best  match 
between  the  response  of  the  model  and  the  neuron  under  study.  Some  approaches  and  ideas 
of  modeling  with  map  (4)  are  discussed  later  in  this  section. 

Before  considering  possible  models  for  synaptic  coupling  among  map-based  neurons, 
let  us  examine  the  response  of  the  map  to  a  rectangular  pulse  of  external  current  7ext(0- 
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This  would  correspond  to  experiments  with  direct  current  injection  through  a  sharp 
microelectrode.  A  pulse  of  current,  in  the  model,  may  be  presented  by  7wext  =  Iext(nT), 
where  T  is  the  period  of  sampling  (0.5  ms).  The  inputs  fn  and  on  can,  in  general,  be  selected 
as  nonlinear  functions  of  7next  or  be  controlled  dynamically  by  7next  with  a  set  of  difference 
equations. 


3.2  Effect  of  Model  Parameters:  Linear  Dependence  on  Input  Current 

Let  us  first  consider  the  simplest  case  when  fn  and  on  are  linear  functions  of  the  external 
current, 

Pn  =  Pei;X\  CJn  =  <7e/„ext,  (5) 

where  the  constants  /3e  and  oe  are  used  to  adjust  the  balance  between  the  two  input  variables. 
The  responses  generated  by  the  map  for  different  values  of  the  pulse  amplitude  Ia  and 
parameter  fe  are  shown  in  Fig.  3.  The  parameters  of  the  map  are  selected  such  that  the 
neuron  is  in  the  silent  regime:  a  =  4.0  and  o  =  —0.4.  The  waveforms  are  computed  for 
the  case  when  the  value  of  7wext  changes  from  zero  in  the  initial  state  to  7next  =  Ia  at  time 
n  =  3,000.  After  1,000  iterations  (^0.5  s),  the  value  of  7next  returns  back  to  zero.  The  value 
of  g e  is  set  to  one.  The  left  panels  illustrate  changes  in  the  response  as  a  function  of  the 
pulse  amplitude  Ia  with  fixed  /3e.  The  right  panels  illustrate  the  influence  of  the  value  of  /3e 
on  the  response  behavior  when  the  amplitude  of  current  pulse  is  fixed  to  Ia  =  0.45. 

The  dynamical  mechanisms  behind  the  neuron  model’s  responses  shown  in  Fig.  3  are 
clear  from  the  phase  portraits  presented  in  Fig.  4,  where  the  curves  Ss  and  Su  correspond  to 
the  stable  and  unstable  branches  of  the  slow  manifolds  in  the  limit  /z  — ►  0.  The  fixed  point, 
O ,  of  the  map  (4)  is  located  at  the  intersection  of  these  curves  with  the  y-equilibrium  line 
x  =  a  —  1 .  When  on  =  0,  the  system  stays  in  the  stable  fixed  point  O  shown  by  the  open 
circle  on  the  stable  branch  Ss.  When  a  depolarizing  pulse  is  applied,  the  y-equilibrium  line 
shifts  up  by  the  distance  D  =  crn  and,  as  a  result,  the  location  of  O  on  the  slow  manifold 


Fig.  3  Waveforms  of  xn  for  the 
responses  of  map  (4)  to  a 
rectangular  pulse  of  amplitude  Ia 
and  duration  1 ,000  iterations 
computed  for  different  pulse 
amplitudes  {left)  and  the  input 
balance  parameter  /3e  {right).  The 
shapes  of  the  pulse  are  shown 
below  the  traces  of  xn.  The 
parameter  values  are  a  =4.0, 
a  =  —0.4,  cre  =  1.0,  and 
a  -  Ia  =  0.630,  /3e  =  0.224; 
b- /a  =  0.517,  /3e  =  0.224; 
c  -  Ia  =  0.439,  /3e  =  0.224; 
d  -Ia  =  0.394,  /3e  =  0.224; 
e  -  =  0.3,  Ia  =  0.45; 

f-  pe  =  0.2,  Ia  =  0.45; 
g-Pe  =  0.1,  Ia  =0.45; 
h  -  pe  =  0.0,  Ia  =  0.45 
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Fig.  4  Phase  portraits  of  the  map  (4)  in  the  plane  of  variables  un  =  yn  +  and  xn  explaining  the  dynamical 
mechanisms  behind  the  main  properties  of  waveforms  for  the  responses  shown  in  Fig.  3.  The  samples 
generated  by  the  map-based  model  are  shown  with  red  dots.  To  follow  the  trajectory  of  the  map,  these 
dots  are  connected  with  straight  yellow  lines,  a,  b,  and  c  Phase  portraits  for  the  trajectories  shown  in  Fig.  3h, 
f,  and  d,  respectively 


changes.  The  location  of  O  during  a  pulse  of  current  is  shown  in  Fig.  4  by  the  black  filled 
circle.  When  /3e  =  0,  the  trajectory  of  the  system  moves  from  the  initial  equilibrium  state 
towards  the  new  one  along  the  manifold.  If  the  perturbed  fixed  point  in  the  new  state  is 
unstable  (i.e.,  located  in  the  branch  Su),  then  the  map  generates  spikes;  see  Fig.  4a.  After  the 
pulse  ends,  the  system  returns  to  the  original  state.  The  corresponding  waveform  is  shown 
in  Fig.  3h. 

The  presence  of  the  input  function  /3n  in  the  first  equation  (4)  dramatically  changes  the 
transient  dynamics  of  the  map;  see  Fig.  4b  and  c.  This  effect  is  easy  to  analyze  using  the 
phase  plane  ( un ,  xn)  instead  of  (yn,  xn).  In  the  case  of  nonzero  /3e ,  the  pulse  of  current  makes 
the  trajectory  jump  parallel  to  axis  un  over  a  distance  /3eIa.  As  a  result,  after  one  iteration,  the 
trajectory  gets  far  into  the  spiking  regime,  following  which  the  system  slowly  approaches 
the  new  state,  a  periodic  motion,  in  the  case  shown  in  Fig.  4b,  or  a  stable  equilibrium, 
Fig.  4c.  When  the  pulse  of  current  turns  off,  the  trajectory  jumps  back  parallel  to  axis  un. 
It  may  overshoot  the  original  state  (resting  potential)  and,  as  a  result,  some  time  is  needed 
before  the  system  returns  to  it,  sliding  along  the  stable  branch  Ss. 

One  of  the  characteristics  of  the  neuronal  response  patterns  that  can  be  tested  in  bio¬ 
logical  experiments  is  the  time-dependence  of  the  frequency  of  spiking  and  the  variation 


45  Springer 


Oscillations  and  synchrony  in  cortical  network  models 


287 


Fig.  5  The  plots  of  fs  vs  n 
computed  for  the  waveforms 
shown  in  Fig.  3a,  b,  c,  and  d. 


'  3000  3200  3400  3600  3800  4000 


of  this  dependence  with  the  amplitude  of  the  depolarizing  pulse  [31].  This  dependence, 
computed  for  the  waveforms  shown  in  the  left  panels  of  Fig.  3,  is  presented  in  Fig.  5. 

The  role  of  the  parameter  /3e  may  be  seen  from  the  traces  obtained  for  a  fixed  value  of 
Ia  =  0.45  and  six  different  values  of  /3e  (see  Fig.  6).  This  parameter  accelerates  the  firing 
at  the  beginning  of  the  injected  pulse  followed  by  a  decrease  in  spiking  frequency  down  to 
the  rate  set  by  the  input  current  on  =  Ia.  The  rate  of  deceleration  is  controlled  by  the  value 
of  the  parameter  /z. 

Following  the  depolarizing  current  pulse,  f$n  changes  from  a  positive  value,  /3eIa ,  to 
zero,  resulting  in  a  sharp  decrease  in  the  value  xn,  known  as  after-hyperpolarization  (see 
Fig.  3e,  f).  For  a  continuously  spiking  neuron,  this  effect  can  appear  as  the  formation  of  an 
interval  of  silence  after  the  depolarizing  pulse  ends  (see  Fig.  7a).  When  a  hyperpolarizing 
pulse  of  current  /next  is  applied  to  a  continuously  spiking  neuron,  the  parameter  /3e  helps  to 
model  the  following  dynamics:  The  pulse  shuts  off  spiking.  When  the  pulse  is  switched  off, 
the  neuron  immediately  starts  firing  with  a  high  frequency.  The  rate  of  spiking  decreases, 
eventually  converging  back  to  normal  (see  Fig.  7b). 

3.3  Effect  of  Model  Parameters:  Nonlinear  Dependence  on  the  Input  Current 

In  order  to  obtain  a  more  realistic  response  behavior  for  modeling  a  particular  type  of 
neuron,  it  can  be  useful  to  consider  a  nonlinear  relation  between  input  functions,  on ,  f}n, 
and  the  injected  current  7^ext.  This  can  be  especially  helpful  for  capturing  the  effects  of 
after-hyperpolarization  and  rebound  bursts.  The  rest  of  this  section  illustrates  briefly  two 
approaches  to  the  design  of  a  dynamical  dependence  of  the  input  functions  on  the  current. 

Consider  a  neuron  demonstrating  the  following  type  of  behavior:  At  baseline,  given  by 
the  resting  potential,  the  neuron  does  not  fire  spikes.  It  starts  to  spike  tonically  when  a 
rectangular  pulse  of  depolarizing  current  is  applied  over  a  limited  duration.  When  the  neuron 


Fig.  6  The  dependence  of  fs  on 
time  n  plotted  for  different  values 
of  fie.  The  other  parameters  of 
the  map  are  the  same  as  in  Fig.  3, 
right  panels 
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Fig.  7  The  map  responses  to  a  depolarizing  (a)  and  hyperpolarizing  (b)  rectangular  pulse  with  a  duration 
of  1,000  iterations  in  the  case  of  a  continuously  spiking  map.  The  parameters  of  the  maps  are  a  =  4.0, 
or  =  4Me  -  2,  and  /3e  =  0.3.  a  Ia  =  0.2  and  b  Ia  =  -0.2 


is  driven  by  a  slow  ramp  of  depolarizing  current,  the  threshold  of  excitation  depends  on  the 
slew  rate.  When  the  slew  rate  is  too  slow,  the  neuron  does  not  spike,  even  at  levels  of  injected 
current  far  beyond  the  expected  values  for  the  spiking  activity  threshold.  Such  behavior, 
known  as  adaptation,  is  exhibited  by  the  GM  neuron  of  the  stomatogastric  ganglion  of  the 
Californian  lobster  [30].  It  is  clear  from  Fig.  1  that,  if  the  value  of  an  is  proportional  to  the 
current,  the  map  cannot  reproduce  such  behavior.  For  values  of  on  where  the  map  crosses 
the  threshold  level  ot h,  it  starts  spiking  for  arbitrarily  low  slew  rates  of  In  order  to  design 
a  map  that  correctly  models  such  a  neuronal  behavior,  one  can  define  the  input  function  on 
with  the  following  equation: 


&n+ 1  —  (1  +  In+ 1 

where  a  small  parameter,  /xa  (0  <  /xa  <<  /x  <<  1),  defines  the  characteristic  rate  of 
adaptation  of  the  map  to  slow  input  an.  In  this  case,  the  maximal  value  of  on  equals 
(In+ 1  —  In)/ i±a-  Therefore,  if  the  parameters  of  the  map  are  in  the  regime  where  the  neuron 
is  silent,  at  finite  distance  from  <jth5  then  there  exists  a  low  bound  for  the  slew  rate  of  the 
ramp  of  current  that  results  in  spiking. 

The  input  function  /3n  may  be  described  by  a  simple  dynamical  model  used  to  adjust  the 
dependence  of  the  spike  frequency  fs  on  the  time  after  onset  of  the  injected  current  pulse. 
Consider  a  simple  example  when 

fin+ 1  —  (1  M /0  fin  T"  In+ 1? 

where  parameter  /x^  defines  a  characteristic  time  delay  to  the  fast  reaction  of  the  map  to  the 
input.  One  can  see  that,  when  a  rectangular  pulse  of  current  is  applied,  /3n  does  not  change 
sharply  but  drifts  towards  the  value  /3eIn+\  exponentially.  As  a  result,  the  beginning  of 
the  first  spike  is  delayed  and  the  spiking  can  first  accelerate  and  then  decelerate.  Numerical 
simulations  indicate  that  this  approach  is  effective  only  when  /x^  is  of  the  order  of  /x.  Notice 
that,  when  /x^  =  1,  the  equation  for  /3n  transforms  into  the  simple  linear  relation  /3n  =  /3eIn. 

After  the  desired  response  of  the  map  to  an  external  current  is  defined,  one  may  use  well- 
known  models  for  different  types  of  synaptic  currents  to  build  a  model  describing  a  network 
of  coupled  neurons. 


4  Modeling  of  Synaptic  Inputs 

With  map-based  models,  the  equation  for  synaptic  currents  can  be  adopted  from  known 
kinetic  models,  but  they  need  to  be  rewritten  in  the  form  of  difference  equations.  Examples 
of  such  equations  were  given  in  [23].  Here,  we  rewrite  these  equation  in  terms  of  synaptic 
conductances. 
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The  simplest  model  of  synaptic  interaction  includes  a  step-like  change  in  the  conductance 
g  of  the  postsynaptic  neuron  by  some  fixed  amount  gsyn  (maximal  conductance)  immedi¬ 
ately  after  a  presynaptic  spike  occurs  followed  by  its  exponential  delay  with  a  time  constant 
r .  An  equivalent  map-based  equation  for  the  conductance  can  be  written  as 


Ci  =  xsT  + 


gsyn,  spikepre, 

0,  otherwise, 


and  the  synaptic  current  is  computed  as 

In  =  ~gn(x%0St  ~  Xrp), 


(6) 

(7) 


where  gsyn  is  the  strength  of  the  synaptic  coupling  and  the  indices  pre  and  post  stand  for 
the  presynaptic  and  postsynaptic  variables,  respectively.  The  first  condition,  “spikepre,”  is 
satisfied  when  >  a  +  y^re  +  /3jfe  or  >  0,  i.e.,  when  the  value  xJJre  is  in  the  right¬ 
most  interval  of  the  function  (2).  It  corresponds  to  the  times  when  presynaptic  spikes  are 
generated.  The  parameter  y  =  1/r  in  (6)  controls  the  relaxation  rate  of  the  synaptic  current 
after  the  presynaptic  spike  is  received  (0  <  y  <  1).  The  parameter  defines  the  reversal 
potential  and,  therefore,  the  type  of  synapse:  excitatory  (x^  =  0)  or  inhibitory  (x^  =  —1.1). 

To  include  the  effects  of  short-term  synaptic  depression  or  facilitation,  one  can  modify 
(6)  by  adding  one  more  variable,  dn ,  that  modulates  the  synaptic  strength  based  on  the 
levels  of  presynaptic  activity.  As  a  result,  the  conductance  may  be  written  in  the  form  of  a 
two-dimensional  map, 

r  syn  7  1  _  J  {V&rc  "T  gsyn^4i>  (1  > 

[8n+"  n+l]~[{YgT ,  1- (i-<?)(i -4,)}- 

where  0  <  ri  <  1  set  the  rate  of  the  depression  and  0  <  q  1 
recovery.  The  initial  conditions  for  (8)  are  {gQyn,  ^o}  =  {0,  1}. 
synaptic  current  is  computed  using  equation  (7). 


spikepre, 

otherwise, 


(8) 


controls  the  rate  of  synapse 
As  in  the  previous  case,  the 


5  Map-based  Modeling  of  Specific  Cell  Types 

Different  classes  of  biological  neurons  are  characterized  by  different  intrinsic  firing  patterns 
that  result  from  a  unique  combination  of  intrinsic  currents  specific  to  a  given  neuron  and 
from  unique  distributions  of  these  channels  across  cell  compartments.  A  common  way  to 
uncover  the  electrical  properties  of  individual  neurons  is  to  study  its  responses  to  constant 
external  inputs.  In  this  section,  we  will  discuss  the  modeling  of  response  patterns  of  several 
major  cell  classes  such  as  regular- spiking  (RS),  fast  spiking  (FS),  intrinsically  bursting  (IB) 
and  low-threshold  spiking  (LTS)  neurons. 

5.1  Model  for  RS  Neuron 

In  [23],  the  model  (4)  was  tuned  to  replicate  firing  patterns  of  a  typical  RS  neuron  when 
it  spikes  in  response  to  a  rectangular  depolarizing  pulse  In.  The  input  parameters,  in  this 
case,  were  proportional  to  In  =  /3eIn  and  crn  =  creIn).  The  parameters  of  the  map-based 
model  of  the  RS  neuron  were  set  to  the  following  values:  a  =  3.65,  a  =  0.06,  /x  =  0.0005, 
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Fig.  8  Intrinsic  neuronal  firing  patterns  in  vivo  and  in  the  models,  a  In  vivo  data,  b  Map-based  model. 
Left  panels  show  a  fast  spiking  (FS)  neuron,  middle  panels  show  a  regular  spiking  (RS)  neuron,  and  right 
panel  show  an  intrinsically  bursting  (IB)  neuron.  The  map-based  model  of  the  FS  cell  is  given  by  (9),  (10) 
with  a  =  3.8,  yrs  =  2.9,  /3hp  =  0.5,  yhp  =  0.6,  ghp  =  0.1,  fe  =  0.1,  and  amplitude  of  rectangular  pulse 
A  =  1.6 E~2  (100%).  The  map-based  models  of  the  RS  and  the  IB  cell  are  given  by  (4),  where  a  =  3.65, 
cr  =  0.06,  p  =  0.0005,  <re  =  1.0,  fe  =  0.133,  and  A  =  3.0 FT2  for  the  RS  cell,  and  a  =  4.1,  a  =  -0.036, 
p  =  0.001,  ae  =  1.0,  f>e  =  0.1,  and  A  =  1.0 IF2  for  the  IB  cell.  The  duration  of  the  depolarizing  pulse  in 
the  map-based  simulation  is  870  iterations,  which  is  approximately  430  ms 


=  0.133,  and  oe  =  1.  These  parameter  values  were  obtained  by  tuning  the  model 
to  replicate  the  deceleration  of  spiking  and  after-hyperpolarization  (Fig.  8,  middle)  in  a 
Hodgkin-Huxley  type  model  of  an  RS  neuron.  The  time-rates  of  these  effects  show  a  strong 
dependence  on  the  value  of  the  parameter  /z,  while  their  strength  may  be  tuned  by  setting 
the  values  of  and  oe .  The  value  of  the  parameter  o  is  selected  to  set  a  proper  baseline 
of  the  neuron  and  a  is  tuned  to  match  the  shape  of  the  action  potential  sequence.  This  map- 
based  model  also  replicates  the  resonance  properties  of  the  RS  neuron  observed  in  in  vivo 
recordings  (see  below  and  [25]). 

5.2  Model  for  FS  Interneuron 

Typical  firing  patterns  of  an  FS  neuron  in  response  to  a  rectangular  pulse  do  not  have  the 
spike  deceleration  effect  but  shows  a  noticeable  hyperpolarization  caused  by  each  generated 
spike.  To  capture  this  hyperpolarization  effect,  the  slow  subsystem  in  the  model  (4)  was 
substituted  with  an  equation  for  the  hyperpolarizing  current  7„hp  generated  by  the  action  of 
each  spike  as  follows  [23]: 


Th P  _  v  hp  rhp 

1  Y 


ghp,  if  xn  occurs  at  the  peak  of  a  spike 
0,  otherwise. 


(9) 
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The  parameter  yhp  controls  the  duration,  rhp  ~  (1  —  yhp)_1,  and  the  parameter  ghp  controls 
the  amplitude  of  the  hyperpolarization  current.  The  model  of  the  FS  neuron  in  this  case  may 
be  written  in  the  following  form: 


Xn+l  =  f(Xn,  Xn-i,yIS  +  /JhP/„hp  +  PeIn),  (10) 

where  /s  is  a  constant  defining  the  resting  state  of  the  model  and  /,!'p  is  a  new  slow  variable 
computed  using  (9).  Equations  9  and  10  with  parameter  values  a=3.8,  yrs  =  -2.9,  phP  = 
0.5,  yh p  =  0.6,  gbp  =  0.1,  and  pe  =  0.1  are  used  to  describe  the  dynamics  of  an  intemeuron 
(IN)  (Fig.  8,  left). 


5.3  Model  for  LTS  Intemeuron 

Another  important  class  of  cortical  interneurons  -  LTS  intemeurons  -  is  characterized  by  the 
rebound  response  to  a  hyperpolarizing  current  pulse  and  a  depolarizing  sag — progressive 
repolarization  of  the  cell  membrane  potential  during  hyperpolarization  (see,  for  example, 
Fig.  2  in  [32]).  These  properties  are  similar  to  the  rebound  response  patterns  found  in 
thalamic  relay  and  thalamic  reticular  neurons  and  have  been  shown  to  be  mediated  by  low- 
threshold  Ca2+  and  hyperpolarization-activated  cation  (I/*)  currents  [33]. 

The  map-based  model  is  capable  of  replicating  rebound  burst  generation.  This  is 
illustrated  in  Fig.  9,  showing  the  waveforms  generated  by  the  map  model  (4)  when  it  is 
hyperpolarized  by  a  rectangular  pulse  of  three  different  durations.  An  important  feature 
here  is  that  the  duration  of  the  rebound  burst  increases  with  an  increase  in  the  time  over 
which  the  system  is  hyperpolarized.  The  rebound  burst  may  be  controlled  by  the  parameter 
/3e  and  may  easily  be  understood  by  an  analysis  of  the  transient  behavior  seen  in  the  phase 
portrait  in  the  ( un ,  xn)  plane  (see,  for  example,  similar  phase  portraits  in  Fig.  4).  When  the 
values  of  /3n  change  fast  from  zero  to  a  large  negative  value,  the  trajectory  jumps  to  the 
left  and  can  overshoot  the  new  steady  state.  Over  the  duration  of  the  hyperpolarizing  pulse, 
the  trajectory  slowly  approaches  the  new  steady  state  and  un  and  xn  slowly  grow.  When  the 
hyperpolarization  pulse  is  over  the  value,  jumps  to  the  positive  direction  of  un  and  the 
trajectory  can  end  up  in  the  spiking  area  where  a  burst  of  transient  spikes  can  be  generated. 

In  the  model  (4),  the  rebound  response  and  spike  adaptation  are  controlled  by  the  same 
input  parameter  /3n  =  /3eIn .  Therefore,  an  increase  in  /3e  will  enhance  both  properties.  To 


Fig.  9  Waveforms  generated  by 
negative  rectangular  current 
pulses  of  amplitude  0.3  and 
durations  of  100,  200,  and  400 
iterations.  The  beginning  and  end 
of  each  pulse  are  shown  by 
arrows.  Parameters  of  the  map 
are  /z  =  0.002,  a  =  3.8, 
a  —  0.15,  /3e  =  0.6,  and 
<xe  =  1.0. 
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control  these  response  features  independently,  we  can  modify  the  relation  between  /3n  and 
In  by  making  it  nonlinear.  For  example, 


bdIn,  if  4  >  0 
j  PkIn,  if  In  <  0. 


(11) 


In  this  model,  spike  adaptation  may  be  controlled  by  the  parameter  /3d  independent  of  the 
rebound  burst,  which  is  controlled  by  fi1' .  Increasing  ft1'  enables  one  to  gradually  enhance 
the  rebound  response  to  a  hyperpolarizing  current  pulse  without  affecting  the  dynamics  of 
the  neuron  in  the  depolarized  state.  Therefore,  in  order  to  build  an  example  of  a  map-based 
model  for  the  LTS  interneuron,  one  may  use  the  model  (4)  with  the  parameters  of  an  RS  cell 
(i.e.,  a  =  3.65,  o  =  0.06,  /z  =  0.0005,  and  oe  =  1)  and  modify  /3n  to  depend  nonlinearly 
on  the  input  current  (11)  with  /3d  =  0.133  and  /3h  =  0.6. 


5.4  Model  for  IB  Cells 


IB  neurons  are  another  important  class  of  cortical  cells.  They  respond  to  depolarizing  inputs 
with  high-frequency  spike  bursts.  In  vivo,  the  intraburst  frequency  is  between  4  and  10  Hz 
[34].  The  number  of  spikes  per  burst  gradually  decreases  over  the  response  duration  [35]. 
In  the  map-based  model  (4),  this  bursting  effect  is  achieved  by  selecting  the  parameter  a 
to  be  larger  than  4,  such  that  the  model  can  produce  bursts  of  spikes;  see  Fig.  1.  Note  that 
the  number  of  spikes  in  each  burst  increases  as  the  values  of  the  parameters  a  and  a  in  the 
parameter  region  corresponding  to  bursting  behavior  increase.  Figure  8,  right  panels,  shows 
the  response  of  the  map-based  model  of  the  IB  cell  obtained  using  (4)  with  the  parameter 
values  a  =  4.1,  o  =  —0.036,  fi  =  0.001,  oe  =  1.0,  and  /3e  =  0.1. 


6  Intrinsic  Resonance  Properties  of  Map-based  Model  Neurons 

Cortical  neurons  display  complex  intrinsic  resonance  properties  that  are  produced  and 
controlled  by  intrinsic  membrane  conductances  that  may  vary  with  the  type  of  neuron. 
During  periodic  stimulation,  neurons  produce  spikes  with  high  reliability  (across  a  set  of 
trials)  for  stimuli  within  a  specific  frequency  range.  The  reliability  of  firing  diminishes  at 
higher  or  lower  frequencies  [36,  37].  To  study  whether  the  map-based  model  is  capable 
of  reproducing  these  important  properties,  sine-wave  modulated  excitatory  synaptic  stimuli 
were  applied  to  the  map-based  model  of  an  RS  neuron  [25].  The  amplitude  of  the  periodic 
input  modulation  was  selected  such  that  a  neuron  remained  silent  for  a  constant  input  equal 
to  the  maximum  of  the  sine-wave  amplitude.  For  low  modulation  frequencies,  both  in 
vitro  experiments  and  the  modeled  neurons  produced  spikes  at  only  a  few  input  maxima, 
making  the  reproducibility  of  spike  trains  from  one  trial  to  another  low  (see  top  panels  in 
Fig.  10).  As  the  input  modulation  frequency  increased,  the  reliability  of  spiking  improved. 
In  a  range  of  moderate  frequencies  (3-13  Hz),  both  in  vitro  and  modeled  neurons  displayed 
almost  perfect  one-to-one  locking  with  the  input  modulation.  The  origin  of  this  resonance 
is  subthreshold  oscillations  observed  in  real  neurons  [36]  and  captured  by  the  dynamics  of 
the  map-based  model  (see  [26]  for  details).  For  very  high  frequencies,  the  reproducibility  of 
spike  trains  was  reduced  again,  reflecting  the  inability  of  both  in  vitro  and  modeled  neurons 
to  follow  high-frequency  inputs  (see  bottom  panel  in  Fig.  10).  Increase  of  a  in  (4)  produced 
an  effect  similar  to  carbachol  application  in  vitro — a  shift  of  the  resonance  peak  to  a  higher 
frequency  [25]. 
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Fig.  10  Spiking  of  RS  neuron 
in  response  to  100  Poisson 
distributed  independent  input 
spike  trains  with  sine -wave 
modulated  mean  rate: 

R  =  fc(l  —  sin (2 fmt )),  where 
fc  =  200  Hz,  fm  is  the  frequency 
of  periodic  modulation. 

Sampling  rate:  one  iteration  = 

0.5  ms.  Six  different  modulation 
frequencies  fm  and  16  trials  for 
each  frequency  (raster)  are 
shown.  Each  dot  represents  a 
spike.  Map-based  model  of  RS 
cell  (left)  and  experiment  with  a 
RS  cell  of  rat  prefrontal  cortex 
(right) 
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7  Oscillations  in  Cortical  Networks 

Fast  network  oscillations  in  the  30-80-Hz  range  are  associated  with  attentiveness  [38] 
and  sensory  perception  [39]  and  have  a  strong  relation  to  both  cognitive  processing  and 
temporal  binding  of  sensory  stimuli  [40,  41].  These  oscillations  are  found  in  different 
brain  systems,  including  the  cerebral  cortex,  hippocampus,  and  olfactory  bulb.  Cortical 
gamma  oscillations  may  become  synchronized  within  1-2  ms  over  distances  of  up  to 
a  few  millimeters  [39]  and  have  been  shown  to  be  critical  for  feature  integration  in 
perception  and  may  control  the  occurrence  and  polarity  of  Hebbian  modifications.  Despite 
the  importance  of  these  processes,  the  genesis  and  synchronization  mechanisms  of  fast  beta 
(15-30  Hz)  and  gamma  (30-80  Hz)  activities  remain  poorly  understood.  In  this  study,  we 
used  computational  network  models  to  study  basic  synaptic  mechanisms  and  properties  of 
fast  network  oscillations. 

7.1  Gamma  Oscillations  and  Synchronization  in  Small  Neuronal  Circuits 

We  first  considered  a  minimal  network  model  consisting  of  three  excitatory  (RS  type)  and 
one  inhibitory  (FS  type)  neuron  (Fig.  11,  left).  The  RS  neurons  produced  tonic  spiking 
with  random  initial  phases  in  the  absence  of  synaptic  coupling.  Neuronal  dynamics  was 
simulated  using  map-based  models  implemented  with  difference  equations  (maps)  for  RS 
and  FS  cells  described  in  Section  5.  The  synaptic  connections  were  modeled  with  equations 
(6)  and  (7),  where  the  strength  of  the  individual  synapses  gsyn  was  scaled  by  the  number 
of  the  presynaptic  cells.  The  parameter  y  for  the  excitatory  and  inhibitory  synapses  was 
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Fig.  11  Oscillations  in  the  model  of  three  PY  neurons  and  one  IN.  Left,  excitatory  (RS-type)  PY  neurons  are 
mutually  connected  with  inhibitory  (FS-type)  IN.  Right,  increase  of  excitatory  input  from  PY  neurons  to  IN 
(from  top  to  bottom )  changed  the  type  of  oscillations.  Different  colors  indicate  three  different  PY  neurons 


set  at  y  =  0.4  and  y  =  0.3,  respectively.  When  excitatory  (AMPA-type)  connections  from 
pyramidal  (PY)  neurons  to  the  IN  and  inhibitory  (GABAa- type)  connections  from  IN  to  PY 
neurons  were  introduced,  the  network  oscillations  became  synchronized  in  the  20-40-Hz 
frequency  range  (Fig.  11,  right).  The  strength  of  the  coupling  between  IN  and  PY  neurons 
determined  the  synchrony  and  the  resonance  mode  of  synchronization.  For  weak  PY-IN 
coupling,  all  three  PY  neurons  fired  together  with  minimal  jitter.  The  IN  typically  fired  a 
few  milliseconds  later.  Increasing  the  strength  of  the  excitatory  input  from  PY  neurons  to 
IN  changed  the  frequency  (/)  of  oscillations  and  also  transformed  oscillations  from  1/1 
resonance  mode  (fpy  =  /in)  to  1/2  mode  (2  fPY  =  /in)  and  then  to  1/3  mode  (3  fpy  =  /in) 
(Fig.  11). 

7.2  Synchronization  Regimes  in  Large  One-dimensional  Network  Models 

We  constructed  a  large-scale  one-dimensional  model  including  layers  of  excitatory  PY  and 
inhibitory  IN  cells  (Fig.  12). 

Mean  field  potential  (FP)  was  estimated  as  the  average  of  membrane  potentials  of  all 
PY  neurons  in  the  network,  and  the  network  activity  was  characterized  by  the  frequency 
of  the  highest  peak  in  the  power  spectrum  of  the  mean  field  oscillations.  Results  of  this 


•  •  • 


Fig.  12  Sketch  of  a  two-layer  network  containing  RS  PY  cells  and  FS  INs.  The  radius  of  connection  footprint 
was  eight  neurons  (16  presynaptic  neurons)  for  AMPA-type  excitatory  PY-IN  synapses  and  two  neurons  (five 
presynaptic  neurons)  for  GABAa- mediated  IN-PY  synapses 
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Fig.  13  Resonance  phase 
locking  in  one-dimensional 
network  model  of  256  PY 
neurons  and  64  IN  neurons. 
Frequency  corresponding  to  the 
highest  peak  in  the  power 
spectrum  of  the  mean  field  of  PY 
cells’  activity  is  plotted  as  a 
function  of  synaptic  coupling 
between  PY  cells  and  INs.  Note 
several  bands  corresponding  to 
different  resonance-locking 
modes  between  FP  and  PY 
neurons.  The  ratio  /py//fp 
changed  from  one  ( lowest  band) 
to  four  {highest  band) 


analysis  are  presented  in  Figs.  13  and  14.  For  a  large  range  of  synaptic  coupling  coefficients, 
the  network  displayed  fast  oscillations  in  the  gamma  (Fig.  13,  left  half-plane)  or  beta 
(Fig.  13,  upper-right  comer)  frequency  ranges.  The  specific  dynamical  state  of  the  network 
varied  significantly  across  the  plane  of  parameters  shown  in  Fig.  13.  Consider  the  ratio 
C/pyA/fp)  of  the  frequency  of  a  single  PY  cell  oscillation  (/py)  to  the  mean  field  frequency 
(/fp).  Variations  of  power  spectra  corresponding  to  an  increase  in  PY-IN  coupling  strength 
presented  in  Fig.  14  shows  that  the  different  phase-locking  modes  were  altered,  changing 
the  frequency  ratio  in  steps.  This  behavior  is  consistent  with  our  study  of  a  small  neuronal 
circuit  (Fig.  11).  Increase  of  PY-IN  coupling  was  followed  by  a  step-like  increase  in  the 
mean  field  frequency  (/?p)  while  the  frequency  remained  almost  unchanged  between  the 
bifurcation  points  (Fig.  14,  top/left).  Also,  the  frequency  of  the  main  harmonic  of  the  PY 
oscillations  (/>y)  decreased  in  the  same  step-like  matter  (Fig.  14,  middle/left).  As  a  result, 
the  frequency  ratio  />yA/fp  changed  at  each  bifurcation  point  (1/1)  —>  (1/2)  —>  (1/3).  The 
frequency  of  the  IN’s  oscillations  was  always  equal  to  the  frequency  of  mean  field  activity 
(Fig.  14,  bottom/left).  When  IN-PY  coupling  was  increased,  this  bifurcation  stmcture 
disappeared  (Fig.  14,  right).  These  results  suggest  that  (1)  input  from  GABAergic  INs  is 
critical  for  synchronizing  the  network  oscillations  in  the  gamma-beta  frequency  range  and 
(2)  the  balance  between  excitatory  (PY-IN)  and  inhibitory  (IN-PY)  coupling  controls  the 
frequency  of  oscillations  and  the  resonance  mode  (/py/ /fp  ratio). 

7.3  Spatiotemporal  Properties  of  Synchronization 

Our  study  predicts  that  oscillations  could  be  globally  synchronized  among  neurons  with 
different  resting  potentials  and  connected  with  synaptic  delays  even  when  each  neuron 
starts  its  evolution  from  a  different  initial  condition.  An  example  of  this  regime  is  shown  in 
Fig.  15a.  Oscillations  in  different  spatial  sites  were  synchronized  within  a  1-2-ms  range 
while  the  direct  propagation  of  activity  between  remote  sites  would  require  more  than 
20  ms.  During  fast  globally  synchronized  oscillations,  the  synchrony  between  PY  neurons 
was  mediated  by  a  very  fast  spread  of  excitation  across  INs.  Excitatory  PY-PY  connections 
were  excluded  from  these  simulations.  The  boundary  between  already  active  and  still  silent 
IN  neurons  propagated  much  faster  than  synaptic  delays  would  allow,  indicating  that  IN- 
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Fig.  14  Spectrograms  of  oscillation  in  one-dimensional  IN-PY  network.  Increase  of  PY-IN  coupling  was 
followed  by  decrease  of  PY  frequency  (/py)  and  increase  of  /fp  and  /in-  Inhibitory  INs  always  fired  at  the 
frequency  of  the  FP  oscillations,  a  giN-PY  =  0.0007.  b  giN-PY  =  0.0015 


PY-IN  synaptic  connectivity  is  not  responsible  for  the  direct  spread  of  activity  from  one 
group  of  neurons  to  another.  Instead,  synaptic  coupling  created  and  maintained  a  regime 
of  phase  synchronization  where  nearby  neurons  fire  with  minimal  delays.  The  phase  delay 
was  small  even  between  remote  sites.  Changing  synaptic  coupling  transformed  the  globally 
synchronized  state  either  to  a  regime  of  propagating  gamma  waves  (Fig.  15b)  or  to  a  regime 
of  asynchronous  oscillations  (Fig.  15c).  The  representation  of  visual  objects  in  the  primary 
visual  cortex  of  awake  monkeys  has  been  found  to  correspond  with  waves  of  gamma  activity 
propagating  across  the  cortical  surface  [42].  It  was  proposed  that  the  phase  continuity  of 
such  gamma- waves  may  be  the  basis  of  spatial  feature  binding  across  entire  objects  [43]. 
Our  study  suggests  that  the  onset  of  global  synchronization  is  possible  when  the  network  is 
homogeneous  enough — all  the  synaptic  weights  and  intrinsic  parameters  of  the  individual 
units  are  selected  within  the  same  region  of  synchronization  (see  (Fig.  13)).  However,  when 
different  units  belong  to  different  synchronization  bands  (e.g.,  when  variability  between 
units  is  large),  the  network  dynamics  leads  to  regimes  with  activity  propagating  from  one 
part  of  the  network  to  another  [e.g.,  in  the  form  of  rotating  spiral  waves  (Fig.  15b)].  This 
creates  a  finite  time  lag  between  oscillations  in  remote  sites  (Fig.  15b,  bottom). 

The  study  suggests  that  FS  inhibitory  INs  mediate  synchronization  of  cortical  exci¬ 
tatory  neurons,  creating  large  areas  of  stable  synchronized  network  oscillations  where 
excitatory  cells  are  phase-locked  to  the  field  at  rational  frequencies.  Changing  the  strength 
of  inhibitory  coupling  controls  transitions  between  persistent  (/py//fp  <<  1)  [44]  and 
transient  (/pyA/fp  ~  1)  [45]  types  of  gamma  oscillations.  Long-range  excitatory  connec¬ 
tions  between  PY  cells  are  not  required  for  global  synchrony.  The  strength  of  synaptic 
interactions  between  PY  cells  and  INs  defines  the  synchronization  state  of  the  network: 
either  global  synchronization  with  near  zero  phase  lag  between  remote  sites  (see,  e.g.,  [39]) 
or  local  synchrony  mediated  by  waves  of  gamma  activity  propagating  through  the  network 
(see,  for  example,  [42]). 
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Fig.  15  Oscillations  in  two-dimensional  network  model  of  256  x  256  PY  neurons  and  128  x  128  INs. 
Snapshots  of  activity  in  inhibitory,  IN,  population  {top)  and  cross  correlation  of  local  FPs  (average 
activity  of  100  x  100  PY  neurons)  between  two  remote  spatial  areas  {bottom)  are  shown,  a  Nearly  global 
synchronization-phase  locking  with  zero  phase  shift,  gpy-iN=0-002,  giN-PY=0.0007.  b  Local  synchro¬ 
nization  mediated  by  two  rotation  spiral  waves,  gpy-iN=0.0015,  gm-?Y =0.0007.  c  Asynchronous  state, 
gPY_iN=0.002,  giN-PY=0.0015.  In  top  panels ,  red ,  depolarizing  (spikes),  and  blue ,  hyperpolarizing  potentials 


8  Conclusion 

In  this  manuscript,  we  presented  a  complete  description  of  discrete  time  neuronal  models 
capable  of  simulating  realistic  firing  patterns  of  many  different  classes  of  biological  neurons. 
Specifically,  response  patterns  of  RS,  IB,  FS,  and  LTS  cortical  cells  and  INs,  as  well 
as  thalamic  relay  neurons,  have  been  simulated.  For  the  class  of  RS  neurons,  we  also 
showed  that  complex  resonance  properties  of  cortical  RS  neurons  during  periodic  sinusoidal 
external  stimulation  can  be  correctly  captured.  Based  on  the  reduced  model  design,  we 
studied  a  large-scale  cortical  network  model  showing  a  spectrum  of  oscillatory  states 
including  fast  gamma  range  (30-80  Hz)  oscillations.  Our  approach,  combining  biological 
feasibility  with  high  computational  efficiency,  may  be  applied  to  study  many  biological 
problems  where  the  size  of  the  network  plays  an  important  role.  Indeed,  there  are  already 
a  number  of  studies  of  neuronal  dynamics  based  on  this  model  [46-51].  The  simple 
mathematical  implementation  of  these  models  makes  them  a  useful  tool  for  education. 
These  models  can  be  used  to  illustrate  the  basics  of  dynamical  mechanisms  behind  specific 
forms  of  firing  patterns  and  network  behavior.  Finally,  the  map-based  models  create  new 
opportunities  for  the  design  of  real-time  systems  that  reproduce  the  functionality  of  real 
biological  networks  containing  a  large  number  of  neurons  that  can  be  used  in  biomimetic 
robots  and  other  applications. 
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